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1. Introduction 

In quantum field theory the perturbative calculation of a given scattering amplitude or cross 
section requires the evaluation of Feynman diagrams. Especially at higher orders this is a difficult 
task and over the last years, a large variety of methods has been devised to deal with the problem of 
calculating Feynman integrals, see e.g. [1] and references therein. The necessary integration over 
the propagators of the virtual or unobserved particles are typically carried out in momentum space 
and divergent integrals are regularized dimensionally by shifting from 4 to D = 4 — 2s dimensions 
of space-time. For a given Feynman integral the main task is then the derivation of an analytical ex- 
pression in terms of known functions with well-defined properties, which at the same time permits 
a Laurent expansion in the small parameter e = (D — 4)/2. 

In these proceedings we want to focus on particular progress in this direction through the 
systematic and efficient approach to solve difference equations. To that end, we start by briefly 
stating the physics case and present the necessary mathematical definitions. Then we provide 
explicit examples from recent higher order perturbative calculations in Quantum Chromodynamics 
(QCD) and finish with a summary and an outlook. 

1.1 Setting the stage 

For a given scattering process the Feynman integrals are classified by the number n of external 
legs (n-point functions), by the number of independent loops / and by the topology of the associated 
graph. Initially, the Feynman integrals will appear as tensor integrals /Mij/fav- where /I,- denote 
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Lorentz indices. Subsequently tensor integrals can be mapped to scalar integrals by numerous 
methods, thus we are dealing with expressions 

/(D;v„...,v„, = 0-« 

where the momenta pu i = I + 1, • • • ,n, are related by energy-momentum conservation pt = 
f(pi,...,pi), and nit denote the masses of the associated particles. The powers of the propa- 
gators are v, and D is the (complex) space-time dimension. In the most general case Feynman 
integrals such as in Eq. ( |1.1| ) may depend on multiple scales, for instance the masses m, but also 
the non-vanishing scalar products of external momenta. 

In the following, we will focus on perturbative calculations to higher orders in massless QCD 
for single-particle inclusive observables. These depend on a single scaling variable x only, with 
x G [0, 1]. Prominent examples are structure functions in deep-inelastic scattering (DIS), fragmen- 
tation functions in electron-positron {e + e~) annihilation or the total cross sections for the Drell-Yan 
process or Higgs production at hadron colliders. All these quantities can be solved directly and sys- 
tematically in Mellin N-space, an approach that has been used in the past within the framework of 
the operator product expansion (and exploiting the optical theorem) in DIS [2,3]. More recently, 
also innovative extensions to other kinematics have been considered [4]. For a generic observable 
G{x), we can write 

l 

dxx N - l ff{x) 

= J dxx N - 1 J dPSW|M({in} — {out})| 2 5(x- /) 

= / dPS«|M({in} - {out})) 2 -^, (1.2) 

where the first line defines the Mellin transform. Subsequently, we express &{x) through the square 
of the scattering amplitude M({in} — ► {out}) for a given set of incoming and outgoing particles 
and dPS*-" 1 ) denotes the integration measure of the Lorentz invariant phase space. The invariant / 
depends on internal and external momenta and introduces the Mellin-Af dependence. As an upshot, 
the observables are mapped to a discrete set of variables (positive integer Mellin N). 







Once the steps in Eq. ( |1.2[ ) are accomplished, the scalar integrals can be reduced algebraically 
to so-called master integrals. The reduction algorithms are based on integration-by-parts, see 
e.g. [1], 

where pi and pj denote any of the loop momenta and the Mellin-A^ dependence is implicit through 
the invariant /, cf. Eq. (|L2|). Upon resolving the constraints from / the reductions in Mellin 
space act on monomials like (pj)~ Vi+N and give rise to systems of linear equations which can be 
solved in terms of a (small) set of master integrals. This step is automatized by using suitable 
(customized) computer algebra programs although in practice limitations arise at higher orders 
through the excessively large size of the systems of linear equations which need to be considered. 
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Due to the explicit dependence on the Mellin variable N the master integrals themselves are 
functions of N. Their functional dependence on ,/V is completely determined by the difference 
equations they satisfy. These difference equations in ,/V are obtained as well from the solutions to 
the integration-by-parts reductions and they are the central topic of the present proceedings. 



With the help of Eq. (1.2) and the solutions to the integration-by -parts reductions, we are thus 



in a position to express a given Feynman integral I(N) in Mellin space in a recursive manner, 

a (N)I(N)+a l (N)I(N-l) + --- + a k (N)I(N-k) =X(N), (1.4) 
which defines a difference equation of order k with the parametric dependence on £ being implicit 



in Eq. ( |l.4[ ). The functions aj(N) are polynomials in ,/V (and £), sometimes they factorize linearly 
in terms of the type N + m + ne with integer m,n. X(N) denotes the inhomogeneous part which 
collects the simpler integrals and in a Laurent series in e it is typically composed of harmonic sums 

Sff, (in = m\,... ,mk) of weight w = \m\ \-\ h \m,]\ and, possibly, in combination with values of 

the £ -function and powers of ,/V. 

From a practical view-point in quantum field theory, an approach like Eq. ( |1.4[ ) allows for easy 
checking at fixed values of the Mellin moment N. The analytical solution to Eq. ( |l~4| ) requires 
concepts and algorithms of symbolic summation (discussed e.g. in [5]) and has been possible in all 
cases we have encountered in QCD calculations [6-9]. However, the recurrence relations can also 
be embedded in a mathematical framework, to which we will turn in the following. 

2. Solving recurrence relations in difference fields 

The Mellin space approach to Feynman integrals at higher orders which we have briefly 
sketched above leads us to the following key problem: 

Given sequences ao(k),... ,a m (k) and fi(k),... ,fd(k), find all constants c\,. .. ,Cd, free of k, 
and all g(k) such that the parameterized linear difference equation 

a m (k)g{k+m)+a m ^\(k)g(k+m-\)-\ \-a (k)g(k) = afi(k)-\ \-c d f d (k) (2.1) 

holds. 

Note that this problem covers various prominent subproblems. Specializing to d = 1, this is 
nothing else than recurrence solving, i.e. we are back to Eq. (|L4|). Taking m = 1 with «o = — 1 and 
a\ = 1 we get parameterized telescoping; in particular, given a bivariate sequence f(n,k), one can 
set fi(k) := f(n + i— l,k) which corresponds in the hypergeometric case to Zeilberger's creative 
telescoping [10]. Finally, choosing d = m = 1 with ao = — 1, a\ = 1, we arrive at telescoping; for 
more details and applications of these summation principles see [11]. 

The algorithmic solution of Eqs. ( |L~4] ) or ( |2.1| ) requires algorithms for symbolic summa- 
tion which are typically implemented in computer algebra systems like Maple, Mathematica or 
Form [12], the latter being most advantageous if large expressions are involved. Specific imple- 
mentations that deal with symbolic summation are e.g. Summer [13] or the recent summation 



package Sigma [14], which can handle problem (2.1), and therefore telescoping, creative telescop 



ing and recurrence solving, if the ao, • • • ,a m and fi,...,fd are given by expressions in terms of 
indefinite nested sums and products. More precisely, Sigma translates such sum-product expres- 



sions into difference fields, the so-called IT£* -fields, and solves the given problem (2.1) there 
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2.1 The construction of difference fields 



Given an equation of the form (2.1), one can represent the coefficients a, and the inhomoge- 



neous parts fi in the so-called IT£* -fields [15]. 

Definition. Let F be a field with characteristic and let a be a field automorphism of F. Then (F, a) 
is called a difference field; the constant field of F is defined by const CT F = {/ G F|cr(/) = /}. A 
difference field (F, a) with constant field K is called a FEZ* -field if F = K(*i ,...,t e ) where for all 
1 < i < e each F, = K(*i, . . . ,?,-) is a transcendental field extension of 1 F,_i = K(*i,. . . ,t,--i) and a 
has the property that <j(ti) = at[ or a(?;) =/, + a for some a G F,-_i. 

Remark. In order to transform, e.g., nested sums in such ITZ* -fields, one exploits the following 
important fact. Suppose we are given a sum T(n) = YJk=\F(k — 1) where we have represented 
F(n) in a ITZ*-field (F,a) with / G F. Then, one can either express T(n) in F by solving the 
telescoping problem: Find g G F with 

o{g)-g = f. (2.2) 

Namely, if we find such a g, we can rephrase T(n) by t := g + c for some c G const CT F, in particular, 
the shift behavior S(n + 1) = S(n) +/(n) is reflected by o(t)=t + f. 

Otherwise, if we fail to find such a g, we adjoin the sum in form of the transcendental field extension 
F(?) where the field automorphism is extended from F to F(f) by the shift behavior a(t) = t + 
/. Then by Karr's remarkable result [15] it follows that the constant field is not enlarged, i.e., 
const ff F(f ) = const CT F. In other words, (F(f), a) is again a IT£*-field. 

In a nutshell, one either can represent a given sum in the already constructed field F by solving 



the telescoping problem (2.2), or otherwise, one can adjoin the sum in form of a transcendental 



extension; for products we refer to [16]. Finally, we emphasize that this construction process is 



completely algorithmically, see Sec. and therefore the translation mechanism can be carried out 



automatically, e.g. in the Mathematica package Sigma. 
2.2 Solving linear difference equations in the ground field 



Suppose we are given Eq. (|2JJ) and we have rephrased the ao,...,a m and f\, . . . ,fd in a IT£*- 
field F = K(t\ ,...,t e );ia short we call such a difference field also the coefficient field of Eq. (|2~l|). 
Then problem ( |2.1| ) can be rephrased as follows: 



Find all c\ , . . . , Cd € K and g G F such that 

a m o{g) + a m _ i a m ~ 1 (g) + • • • + a g = c x f\ + • • • + c d f d . (2.3) 

Remark. Eq. ( |2.3| ) can be solved in F by solving several such problems in the subfield G := 
K(fi,. . . Namely, we arrive at the following reduction [17]; we set t := t e , i.e., F = G(t): 

Reduction I (denominator bounding). Compute a nonzero polynomial d G G[t] such that for all 



ci G K and g G G(t) with Eq. (2.3) we have dg G G[t]. Then it follows that 



?« , + " + 5 %«"<rt =«,* + ... + „/, (2-4) 



'We set Fi 



; 
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for g' G G[t] if and only if Eq. (Q with g = g'/d. 

Reduction II (degree bounding). Given such a denominator bound, it suffices to look only for 
c\ G K and polynomial solutions g G G[t] with Eq. (|23|). Next, we compute a degree bound b G No 
for these polynomial solutions 

Reduction III (polynomial degree reduction). Given such a degree bound one looks for a G K and 
gi G G such that Eq. Q holds for g = £f =0 g/. This can be achieved as follows. First derive 
the possible leading coefficients gi, by solving a specific parameterized linear difference equation 



given in G, then plug its solutions into Eq. (2.3) and recursively look for the remaining solutions 
g = Yd=o Sit'- Thus one can derive the solutions of Eq. (2.3) given in G(t) by solving several such 



parameterized linear difference equations given in the ground field G. 

Together with results from [15, 18-20] it has been shown in [17] that this reduction leads to 



a complete algorithm for problem ( |2.3| ) with m = 1 ; as result we obtain a streamlined and simpli- 



fied version of Karr's algorithm [15] whose main interest is to solve the telescoping problem (|2.2|). 
Moreover, it is shown that this reduction also delivers a method that eventually produces all solu- 
tions for recurrences of higher order m>2. 

2.3 d'Alembertian solution 



Now suppose that we have rephrased Eq. ( |2.1| ) in the form ( [2.3D for a particular coefficient 
field F. Then F is usually too small to find all the required solutions. Therefore, the following 
more general approach is helpful: 
Find all solutions of the form 

h(n) £ b!(h) £ b 2 {k 2 ) ... £ b s {k s ) (2.5) 

,t 1= k 2 =0 k s =0 

where the bi(kj) and h{n) are represented in F or they are defined as products over elements from F. 
Such type of solutions are called d'Alembertian solutions [21], which are a subclass of Liouvillian 
solutions [22]. 

Remark. The d'Alembertian solutions are obtained by factorizing the linear difference equation as 
much as possible into first order linear right factors over the given difference field/ring. Then each 
factor corresponds basically to one indefinite summation quantifier; see [21,23]. We remark that 
finding such a linear right hand factor is equivalent to finding a product solution of the recurrence in 
its coefficient domain. For the rational case, i.e., aj(k),fi(k) G K(k), this problem has been solved 



in [24]. A general version for the IT£* -field situation, which uses the methods described in Sec. 2.2 
as subroutines, has been developed in the summation package Sigma; see also [14,23]. 

We stress the following important aspect: one usually needs simplified representations of the 



solutions ( |2.5| ) for further treatment. If the solutions are given in form of harmonic sums, one can 
derive compact representations by using the algebraic relations of harmonic sums (see e.g. [25]), 
or alternatively, use Sigma for this task. Ongoing research (see e.g. [26]) is dedicated to the 



computation of sum representations of the type (2.5) with optimal nested depth 



It is worth pointing out, that within the formalism of IT£* -fields, a number of exten- 
sions/generalizations can be treated (e.g. in Sigma). These include algebraic objects like (— 1)". 
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Note that such elements cannot be represented in fields, but only in rings with zero-divisors, like 
(1 + (— 1)")(1 — (—1)") = 0; for more details see [23]. One can also consider radical objects like 
y/n or free/generic sequences X n , for which the corresponding algorithms are presented in [27] and 
in [28,29], respectively. 

Summarizing, given an equation of the form (2T ), the mathematical framework of ITZ* -fields 
puts us in a position to solve the corresponding recurrence relations in terms of indefinite nested 
sums and products over the given coefficient field F. This class covers harmonic sums SV, and is 
therefore sufficient for solutions of difference equations for single scale Feynman integrals in the 
Mellin space approach. 



3. Examples of Feynman integrals 

Here we present two examples from QCD calculations to two- and three loops for DIS structure 
functions for single hadron inclusive e + e~ -annihilation [6-9]. The relevant diagrams are displayed 
in Fig. |l| and the respective difference equations for the Feynman integrals I(N) are of second and 
third order. They could be solved by matching to an ansatz of the type 

/(/V) = Y,<j,m)e j S ifl (N)+Y,b(j,fn,k)ejCkSdN), (3.1) 
were the unknown coefficients a and b were then determined with the Summer package [13] by 



inserting Eq. ( J3- lp in the recursion relation ( |1.4[ ). This approach, of course, rests entirely on the 
fact that the solution for the Feynman integrals I(N) as a Laurent series in e is within the space 
spanned by harmonic sums Sg t (N) and Riemann ^-values Q. For all Feynman integrals I(N) which 
were determined by higher order difference equations we found this condition to be fulfilled. 



An alternative path to the solution, based on the mathematical framework developed in Sec. 2.2 
above, is provided by the Sigma package. This we want to discuss next. 

3.1 e + e -annihilation 

The example which stems from single hadron inclusive e + e~ -annihilation is a phase space 
integral for a decay process 1 — > 4 in one-particle inclusive kinematics, see Fig. [j]on the left. The 
scalar diagram describes the decay of a particle with initial momentum q according to q — » p\ + 
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Pi + P3 + P and the dashed vertical line denotes the final-state cut for the production of particles. 
Our example is the master integral R2 (N) from the so-called real-real-emission and it depends on 
the dimensionless variable x E which is the scaled momentum fraction 

x E = 7 ^—t- , < x E < 1 . (3.2) 

q-q 

The solution of the integration-by-parts identities leads to a difference equation of second 
order [9], which reflects the underlying symmetries of the Feynman diagram. We find 



(N + 1 - 2e) (N + 2 - 6e)R 2 (N) - (N - 1 ) {N - 4e)R 2 (N - 2) 
(1 - e)(l - 3e)(2- 3e)(2iV+ 1 - 6e) 
(N-3e)(N+l-3e) 



Ri(tf-2), (3.3) 



where Ri denotes the inhomogeneous part. For a complete solution, we also have to provide the 
initial conditions R2(0) and R2(l)- The explicit expressions are rather lengthy, in particular the 
higher powers in the Laurent series in e and we refer the reader to [9] for details. 

We have provided Sigma with the difference equation (^3|) and the expression for Ri in terms 
of harmonic sums. Subsequently Sigma was able to solve for R2(./V) and we have found complete 
agreement with the published result after insertion of the initial conditions R2(0) and R2(l)- 



3.2 Deep-inelastic scattering 

During the computation of the three-loop QCD corrections to the DIS structure functions 
and Fl we have encountered the Feynman integral shown in Fig. [T] (right). The scalar diagram 
in the DIS case has external momenta p and q and is of the non-planar topology NO22 with unit 
powers of the propagators. Due to the framework of the operator product expansion and the optical 
theorem applied in DIS, we are specifically interested in the imaginary part of N022- In momentum 
space it depends on the scaling variable xg, 

*b = , 0<x B <\. (3.4) 

In Mellin space the associated difference equation is of third order and very difficult [8]. It has 
been obtained by means of a systematic solution of the integration-by-parts identities and is given 
by 

-NO22 (N)N(N - 2) (3N - 5) (N 2 - 1 ) 

+N0 22 (Ar - 1 )N (3N 4 - 2(W 3 + 45N 2 - 40N + 12) 

+N0 22 (Ar - 2)N(3N 4 - UN 3 + 35N 2 - 31N+ 10) 

-N0 22 {N-3){N-l){N-2) 2 (3N-2)(N-3) = 2X(N) , (3.5) 

where we have put e = 0, since the integral NO22 does not contain any divergences. It is completely 
finite and to the accuracy needed in the three-loop calculation [8] positive powers in e were not 
necessary. 
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For the solution of Eq. (3.5), we have to provide the inhomogeneous part X (which is also 
finite), 

v 1 b N-3 N-2 v ' N-3 (N-3) 2 

N-2 (N-2) 2 v ; N-2 (N-2) 2 

- ,2*^- ') + 8«^i> - 12S 2 (») - 16*<^ - 8^i> - «,(») 
N-l (N-l) 2 N-2 N-l V ; 

r \ r i i 1 C 

+84-7^- +24— —-4-^— + 92— — ^ + 56— - — r -48—— +4- 



N-3 N-2 N-2 (N-2) 2 (N-2) 3 N-l N-l 

- m W^~ 36 W^ +4 *N +16 ^ 

S^V-3) ^ S_ 3 (N-2) 5- 3 (iV-l) ,. A 
-56 24 + 8 — 16S_ 3 (N) - 48S_ 3 (N)N 



N-3 (N-3) 2 N-2 (N-2) 2 N-l 

' S " 2 ( A ' ~ 1 ) 165_ 2 (N) - 48S_ 2 (N)N - 84-^- - 24— i— - 36- ^ 



(N-l) 2 v ; v ; N-2 N-2 

-92- — l —— - 56- — 3—5- + 12-2?- + 12- — L_ + 20- — ?— - ? - 12-^ - 20-^ 

(N-2) 2 (N-2) 3 N-l (N-l) 2 (N-l) 3 N 2 N 3 

(3-6) 

along with the initial conditions NO22(0), N022(l) and N022(2). Here, Q denotes the Riemann 
^-function at value i. 

We have again used Sigma for solving the difference equation (3.5), given the expression for 
X from Eq. (3^6) and Sigma has provided us with the correct analytical answer for NO22 once also 
the initial conditions NO22(0), NC>22(1) and N022(2) had been supplied. The solution is given by 
the following very compact expression, 

N ° 22(A ° nWTT) 

8(l + (-l)*)S_ s (tf) 4(l + (-lf)5 5 (N) t ^ N) {4S 3 (N) 4£ 3 



+5-3 (N) 



N+l N+l v ' V N+l N+l 

/ 8 ( 1 + (- If) 4 (2 + (- if) S- 2 (N) 4S 2 (N) \ 
1 N 2 (N+1) N+l N+l I 



/ -12C 3 N 3 + (-lf (8-8N 3 C 3 )+8 ^(-lf^A 
21 ' 1 N 3 (N+1) N+l J 

85_ 3: _ 2 (N) (4-4(-l)")S_ 3 , 2 (AQ (4 + 12(-l) jV )^ 2 , 3 (N) 8%(N) 
AT+1 N+l N+l N+l ' 

(3.7) 

where we have eliminated a number of nested sums in favor of products of harmonic sums through 
their product algebra, leaving only four independent sums of depth two. 
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4. Summary 

In these proceedings, we have presented a short introduction to the problem of calculating 
Feynman integrals at multiple loops. For single-scale problems we have sketched how to formulate 
the problem in Mellin space, how to obtain algebraic reductions for the loop integrals and, finally, 
how these reductions lead to difference equations in Mellin space. While the underlying physics 
problem usually places certain restrictions on the difference equations, which allow for their solu- 
tions to be expressed in terms of harmonic sums, and hence to be constructed with a corresponding 
ansatz, it is advantageous to reconsider the problem in a more general mathematical framework. 

We have reported on progress in this direction by exploiting algebraic properties of differ- 
ence fields, the so-called TIL* -fields. The algorithmic construction of analytical solutions to the 
recursion relations for Feynman integrals over difference fields as, for instance, implemented in the 
Mathematica package Sigma provides great calculational advantages. To that end, we have tested 
the package Sigma on two examples of Feynman integrals which had appeared in recent perturba- 
tive higher-order QCD computations for e + e~ -annihilation and deep-inelastic scattering. With the 
improved mathematical machinery both examples have been re-evaluated with the package Sigma 
and we have found agreement. 

In the future more generalizations or extensions are conceivable. From the need to consider 
physics problems and Feynman integrals depending on multiple scales one arrives at generalized 
sums [30], which in turn are connected to multiple and harmonic poly logarithms [31,32]. A related 
problem is the expansion of (generalized) hypergeometric functions around integer or rational val- 



ues in a small parameter These problems lead to recursion relations similar to Eqs. (1.4) or (2.1) 
although with functional or parametric dependences on more variables. Some algorithms and im- 
plementations exist [33-36] and it will be interesting to pursue this line of research further within 
the framework of IT£* -fields. 
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